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1.  Introduction 

In  this  paper  we  consider  a class  of  methods  which  are  designed 
to  Increase  the  computational  efficiency  of  computer  simulations.  This 
increased  efficiency  is  obtained  by  simulating  a stochastic  process  which 
is  related  to^  but  different  from^  the  stochastic  process  of  Interest. 

Each  event  in  the  new  process  will  correspond  to  several  events  in  the 
original  process.  The  simulation  of  this  new  process  will  then^  in  some 
sense^  proceed  at  a faster  rate  than  that  of  the  original  process. 

The  types  of  processes  for  which  this  technique  may  be  applied  are 
positive  recurrent  Markov  chains  in  both  discrete  and  continuous  time  as 
well  as  semi-Markov  processes.  These  are  all  examples  of  regenerative 
processes  (see  Cinlar  (I975)  or  Crane  and  Iglehart  (1975)),  ^nd  the  method 
to  be  proposed  relies  heavily  on  this  fact.  Under  the  regenerative  assump 
tion  a single  simulation  run  may  be  broken  up  into  randomly  spaced  i.i.d. 
(independent  and  identically  distributed)  blocks^  or  cycles.  This 
allows  the  techniques  of  classical  statistics  to  be  applied  in  analyzing 
the  output  of  the  simulation. 

One  of  the  main  difficulties  with  regenerative  simulations  is  that 
even  though  it  may  be  known  that  the  process  being  simulated  is  regenerati 
the  regenerations  may  be  few  and  far  between.  Thus  even  for  very  long 
simulation  runs^  only  a relatively  few  i.i.d.  cycles  are  observed.  In 
this  case  it  becomes  difficult  to  form  reliable  point  and  Interval 
estimates.  One  possible  remedy  to  this  problem  is  the  use  of  approximate 
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regenerations  (see  Crane  and  Iglehart  (1975a)  or  Gunther  (1975)).  The 
idea  here  is  to  break  the  simulation  run  up  into  blocks  which  are  "almost" 
i.i.d.  by  defining  a sequence  of  "almost"  regenerations.  These  blocks 
are  then  treated  as  if  they  are  i.i.d.  (although  they  are  not).  Since 
"almost"  regenerations  are  defined  so  that  they  occur  more  frequently 
than  actual  regenerations,  one  obtains  more  blocks  using  the  approximation 
than  is  otherwise  possible.  Presumably  this  facilitates  the  formation  of 
point  and  interval  estimates,  although  this  has  never  been  demonstrated 
satisfactorily.  Because  this  method  lacks  a strong  theoretical  founda- 
tion, it  should  be  used  with  caution. 

The  method  proposed  here  also  seeks  to  increase  the  number  of 
blocks  obtained  during  a simulation,  but  does  so  without  resorting  to 
any  approximations.  Instead,  a new  Markov  chain  is  simulated  for  which 
regenerations  occur  more  often  than  for  the  original  chain.  This  new 
chain  is  constructed  In  such  a way  that  point  estimates  and  confidence 
intervals  may  still  be  formed  for  certain  parameters  of  the  original 


chain. 


2.  Derivation  of  the  Modified  Markov  Chain 


We  now  state  this  basic  problem.  Let  {X  , n > O)  be  an  irreducible 

n 

aperiodic,  positive  recurrent  Markov  chain  with  state  space  E = (0,1,2,...) 

and  transition  matrix  P = (p^^  : i,j  £ E).  It  is  then  well  known  that 

there  exists  a probability  distribution  tt  = (tt^  : i £ E)  on  E and  a 

random  variable  X,  having  distribution  tt,  such  that  X^  X (=^  denotes 

weak  convergence),  tt  is  called  the  stationary  distribution  of  the  process 

and  is  usually  unknown  or  difficult  to  calculate.  Let  f be  a real 

valued  function  on  E and  define  r = irf  = E(  f(X)  ] = Y.  TT . f ( i) . We 

i£E  ^ 

shall  be  Interested  in  estimating  r. 

Pick  some  state,  say  0,  in  E and  let  = 0.  Define 


T = inf{n  > T , : X =0) 
m m- 1 n 


m > 1. 


We  say  that  a regeneration  occurs  at 

and  T -1  is  referred  to  as  the  mth 
m 

T is  the  length  of  the  mth  cycle. 

Ttl 


time  T and  the  time  between  T , 
m m- 1 

cycle.  If  T = T - T ,,  then 
■'  m m m-1’ 

Let 


Y 

m 


T -1 
m 


n=T 


S £(x„) 

m-1 


m > 1 . 


Because  (X  , n > 0)  is  positive  recurrent  E,[t  1 = E[t  |X  = i]  < ». 
"■n'”  Im  mu 

If  X = 0 then  {(Y  , T ),  m > 1)  are  i.i.d.  It  is  also  known  that  if 

Trlfl  = Z TT , I f(  i)  I < ®,  then 
i€E 
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.. 


(£.1) 


(see  Crane  and  Iglehart  (1975)).  Let 


Z = Y - rr 
in  m m 


2 2 

By  equation  (2.1).  E-[Z  ] = 0.  Let  a = E-[Z  ] and  assume  that 

Om 


0 < a <00.  Define 


(2.2) 


(2.3) 


^ M M 

r(M)  - Z V \ \ • 

m=l  in=l 


x(N)  = Z f(X  )/N+l 
n=0 


Then  r(M)  -»  r a.s.  (almost  surely)  as  M -»  » and  x(N)  -»  r a.s.  as 
N 00.  Because  {Z  m > 1)  are  i.i.d.  (assuming  X = 0)  it  is  easy 
to  prove  the  central  limit  theorems 


(2.M 


(2.5) 


:x(N)  - r 


a/Eo[Ti] 


172  ^ N(0,1) 


as  M -» 00 


as  N 


where  N(O^l)  is  a normally  distributed  random  variable  with  mean  0 
and  variance  1.  Point  estimates  for  r can  be  given  by  either  (2.2) 
or  (2.3)  and  confidence  intervals  for  r can  be  based  on  (2.4)  or  (2,5) 


t 


We  now  show  how  to  define  a modified  Markov  chain  which  will  have 


t 


a shorter  expected  cycle  length  and  from  which  point  estimates  and 
confidence  intervals  for  r may  be  derived.  Introduce  a new  state 
and  let  = E U A.  The  addition  of  the  state  A will  enable  us  to 
use  Dynkin’s  formula,  which  is  the  basis  of  the  method.  Define  a new 
transition  matrix  P = (p^^j  : i,j  C where 

j = 0 

i € E,  j € E,  j / 0 
i € E,  j = A 
i = j = A 


P is  then  basically  the  same  as  P,  the  major  difference  being  that 
column  0 of  P has  been  placed  in  column  A of  P and  column  0 
of  P is  set  equal  to  0.  If  E = {0,  1,  ...,  m),  then  P may  be 
written  as 

0 1 • • • m A 


Following  the  notation  of  Hordljk,  Iglehart  and  Schassberger  ( I976) 

denote  the  submatrlx  of  P consisting  of  all  rows  and  columns  except 

for  row  and  column  A by  ^P.  Let  n > 0}  be  a Markov  chain 

A 

with  state  space  E and  transition  matrix  P.  For  this  Markov  chain 
each  state  In  E Is  transient  and  A Is  an  absorbing  state.  Extend 
the  definition  of  any  vector  (or  function)  g on  E to  by  setting 

g(A)  = 0.  Define  the  absorption  time  t by 


T = Inf (n  > 0 : = A) 


and  let 


7 - V 

n=0 


Notice  that  If  X^  = X^  a.s.  then  the  pairs  and  (Y^  t) 

have  the  same  distribution  ( In  fact  If  we  simulated  these  two  processes 
using  the  same  stream  of  random  numbers  then 

y(l)  = E^[Yj^]  = E^[Y]  and  t(  1)  = = E^[t].  By  equation  (2.1) 

r = y(0)/t(0).  Since  X^  = A for  all  n > t and  f(A)  = 0 then 


Y = Z f(X„) 
n=0 


By  taking  expectations  we  find  that 


(2.6) 


y = I p"f 
n=0 


j € E 


j , 


and  a^j  is  the  expected  number  of  visits  of  the  Markov  chain  n > O) 

to  state  j given  that  = i.  If  we  interpret  » 0 = 0 then  (2.6) 
and  (2.7)  may  be  written  as  y = Af  and  t = Ae. 

The  derivation  of  the  new  Markov  chain  to  be  simulated  will  be 
based  on  Dynkin's  formula^  which  is  given  in  the  following  proposition. 

For  the  definition  of  a stopping  time  sea  page  118  of  (Jinlar  ( 1975)* 


(2.8)  PROPOSITION.  Let  S be  any  stopping  time  for  the  Markov  chain 


{X  , n > 0)  such  that  d( i)  = E[slx  = i]  < « for  each  i. 
n’  ~ o 


rS-  1 -1 

i)  = E 5:  f(X  )lx  = iU  E[y(Xg)lx  = i] 

Ln=0  J 


If  Tlfl  < 
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Let  fC  n > 0}  be  a Markov  chain  with  state  space  E and  transition 
n’  — 

matrix  R.  This  Markov  chain  also  has  A as  an  absorbing  state.  Let 


T = inf{n  ^ 0 : = A)  ^ 


f - z'  l'(5„)  , 

n=0 


and 


5 . ’i'  d(c„) 

n=0 


Since  E(Y}  = X.  equation  (2.9)  implies  that  E[y|c  = i]  = y( i) . 

n=0 

Similarly  EfSlc^  = i]  = t(i). 

We  now  remove  the  state  A.  Assume  that  = P{Xg  = ® 

for  all  i (this  will  be  the  case  for  all  interesting  stopping  times  S) . 
Define  the  transition  matrix  R = {r^.  ; i^j  £ E}  where 


r.  . 

ij 


j = 0 

j 0 . 


R is  obtained  from  R by  placing  column  A of  R in  column  0 
of  R,  deleting  row  and  column  A of  R'^  and  leaving  R otherwise 
unchanged.  Now  let  fC  , n > 0)  be  a Markov  chain  with  transition 
matrix  R,  Define  = 0^ 


T'  = inf(n  >T':C=0),  m>l, 

m m- 1 n ^ - * 


t’  = T’  - T',  , 
m m m- 1 ' 


m > 1 j 
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T^-1 

Y*  = E h(C  ) , 

n=T-  , 
m-  1 


T'-l 

&•  = "Z  d(Cj  , 


n=T' 


and 


m-  1 


Z • = Y ' - r6  ’ , 
m tn  jn  ' 


m > 1 , 


tn  > 1 , 


m > 1 . 


As  before  if  a.s.,  then  (Y|,  6^,  t|)  has  the  same 

distribution  as  (Y^  5^  t).  In  particular  E[Y|ICq  = i]  = y(i)  = EfYj^lX^  = i] 

and  E[&^(Cq  = i]  = t(i)  = = i],  so  that  E^fZ^l  = 0 Furthermore 

the  times  (T*.  m > 1}  are  regenerations  for  (C  , n > 0}  so  that  if 
tn'  — n 

= 0.  {(Y',  5',  t',  Z'),  m > 1)  are  i.i.d. 

0 » '•'  m’  m»  m’  m''  — 

The  Markov  chain  fC  . n > 0}  is  the  one  that  will  finally  be 

n'  — 

simulated.  The  use  of  the  state  A in  defining  the  stopping  time  S 
automatically  ensures  that  the  return  state  0 cannot  be  inadvertently 
"skipped"  during  the  simulation  of  n > 0} . 

Defining 


(2.10) 

and 

(2.11) 


M M 

^•(M)  = E^  Yy  E 6; 


m=l 


m=l 


, N N 

X (N)  = E h(C  )/  E d(C^)  , 
n=0  n=0 


then  r'(M)  -»r  a.s.  as  M->oo  and  x*(N)->r  a.s.  as  N->'xi,  Let 

,.2 

I ann  >«v>ki  i i ■ m l t i >».  i / 

0 m 

derive  the  central  limit  theorems; 


= E-fZ'^]  and  assume  that  0 < It  is  straightforward  to 
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(2.12) 


£LS  M — > 00  ^ 


Point  estimates  for  r are  now  given  by  either  (2.10)  or  (2.11)  and 
confidence  intervals  can  be  based  on  (2.12)  or  (2.13).  Once  the  transi- 
tion mtrix  R and  functions  h and  d have  been  calculated^  the 
formation  of  point  estimates  and  confidence  intervals  is  essentially  the 
same  as  in  the  regenerative  method. 


3 . Numerical  Discussion 


To  determine  the  amount  of  variance  reduction  obtained  by  simulating 

the  chain  n > 0}  rather  than  n > 0)  we  need  only  compare 

the  variance  terms  in  the  central  limit  theorems  (2.12)  and  (2.5).  This 

2 

variance  reduction^  which  will  be  denoted  by  Rg^  is 


which^  since  Eq[8^]  = simplifies  to 


2 


If  o'  and  a are  approximately  equal  then  Rg  = 

ratio  of  the  expected  cycle  lengths  for  the  two  processes.  This  suggests 
that  S be  made  as  large  as  possible^  however  by  doing  so  the  amount  of 
work  needed  to  compute  R,  h and  d may  be  prohibitive.  Each  transition 
in  the  chain  (C^,  n > 0}  corresponds  to  E[s|c^]  transitions  in  the 
chain  {X  n > 0}  so  that  we  expect  transitions  for  {C  , n > 0)  to 
be  relatively  expensive  to  generate.  The  stopping  time  S should  there- 
fore be  chosen  so  that  R^  h and  d may  be  readily  computed  and  the 
sample  paths  of  n > 0)  may  be  easily  generated.  We  give  several 

examples  of  such  stopping  times. 


(3.1)  EXAMPLE.  Constant  S. 

m- 1 m- 1 

If  S = a constant^  then  h = Z = E 

n=0  n=0 

R is  given  by 


O^ij  > 


‘ ■ .1,  0^-  ’ 


j ^ 0 


j = 0 . 


If  P is  relatively  sparse  the  work  involved  in  computing  h and  d 
should  not  be  too  great  for  small  values  of  m.  If  S = then  R = P^ 
h = f and  d = e,  so  that  this  choice  of  S reduces  to  straighforward 
simulation  of  the  original  Markov  chain  n > 0}.  If  m is  small 

and  if  p^^  z=  0 for  most  i £ E then 


m- 1 m-  1 

d = Z *=  Z P”e  = me 


Since 


E[t^1  = Z pP^e  = Z gR^d  « m Z gR^e  = mEfT’l  , 

n=0  n=0  n=0 


we  expect  that  Rg  « l/ra.  Tables  1 and  2 bear  out  this  speculation  for 
two  birth  and  death  processes^  the  finite  capacity  M/M/ 1 queue  and  the 
repairman  problem.  The  repairman  problem  has  birth  and  death  parameters 


15 


n>\  ^ 0 < i < s 

(n+s-l)X  , s < t < n+s 

^ , 1 < i < c 

= I 

Icn  , c < i < s+n  , 

where  n is  the  number  of  operating  units^  s is  the  number  of  spare 
units,  c is  the  number  of  repairman,  and  \ and  ^ are  the  failure  and 
repair  rates  respectively  of  the  units.  These  continuous  time  problems 
have  been  transformed  into  discrete  time  using  the  methods  of  Hordijk, 
Iglehart  and  Schassberger  ( I976) . 

(3.2)  EXAMPLE.  Exit  Times  of  Sets 

Suppose  E is  partitioned  into  N„  disjoint  blocks  B.,-  i.e.. 

D 1 ' 

Nb 

E = U and  5^08^=  0,  the  empty  set,  for  i ^ j.  Let  |B^! 

denote  the  number  of  elements  in  B.  and  assume  that  Ib  I < ».  Define 
the  stopping  time  Sg  to  be  the  first  exit  time  from  the  initial  block; 
i.e. , 

Sg  = lnf(n  > 0 : 0 B^)  , for  6 B^  . 

Let  f^  = {f(  j)  : j e B.),  e^  = (e(j)  : j € B^}  and  ^P^^  = : k € B., 

£ € ®j}-  denote  the  identity  mar  .ix  of  size  jB^^l  by  jB^j. 

Then 
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and 


-1 


h.  = (I..  - „P.,)  f. 

1 11  0 li'  1 


-1 


= (I, . - ^P. .)■  e. 
i ' ii  0 ii'  1 


O^ij 


i = J 
i j 


(I  - P ) P 
^ ii  0 ii^  0 ij  > 

where  h.  = {h(j)  : j € B.},  d^  = {d( j)  : j G B^)  and  : k £ B^, 

£ € ®j}*  transition  matrix  R is  found  by  setting 


. 

0 ij  * 

j ^ 0 

^ o’^ik  > 

j = 0 

ij 


The  matrix  ^R  is  the  block  Jacob  matrix  for  solving  the  set  of  equations 

y = f + Q^y  (see  Varga  (I962)).  The  Markov  chain  {C^,  n > 0}  is  the 

same  one  that  appears  when  using  the  techniques  of  Heidelberger  ( I978) . 

2 

Tables  1 and  2 give  R„  and  R for  the  finite  capacity  M/m/ 1 queue 

d O 

and  the  repairman  problem  respectively. 


(5.5)  EXAMPLE.  Continuous  Time  Markov  Chains  and  Semi-Markov  Processes. 

By  applying  this  technique  to  continuous  time  Markov  chains^  the  dis- 
crete time  methods  of  Hordijk^  Iglehart  and  Schassberger  ( I976)  are  obtained 
as  a special  case.  Let  {X^,  t > 0)  be  an  irreducible^  positive  recurrent 
continuous  time  Markov  chain  with  infinitesimal  matrix  Q = (l-j  • i,j  £ E). 
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Set  q.  = Z q . . (‘I-  = Dynkin's  formula  remains  valid  for  continuous 

1 ij  1 11 

time  Markov  chains  (with  an  integral  replacing  the  sum  in  (2.18)).  Define 
the  stopping  time  S = inf{t  >0  : ^ ^0'"  then  simulate  a discrete 

time  Markov  chain  n > 0)  with  i j,  h(  i)  = f(i)/q^ 

and  d( i)  = l/q^.  Hordijk,  Iglehart  and  Schassberger  ( I976)  have  shown 
that  in  this  case  a’  < a.  These  techniques  may  also  be  used  to  reduce 
semi-Markov  processes  to  discrete  time.  Once  in  discrete  time,  the  method 
may  be  applied  again  to  obtain  further  variance  reductions. 

!i 


TABLE  1 


V 


f 


i 

► 

& 


2 

Variance  Reductions,  Rg  and  Rg,  for  Finite  Capacity  M/lV  1 Queue 
Obtained  by  Simulating  the  Modified  Markov  Chain: 
r = E(X),  Return  State  = 0 


p 

Stopping  Time 

2 3 

s 

* 

^B 

.25 

.4264 

.2923 

.2942 

.6530 

.5406 

.5424 

.50 

.4590 

.3620 

.4441 

.6775 

.6016 

.666k 

.75 

.4841 

.3585 

.4057 

.6957 

.5987 

.6353 

.90 

.4883 

.3422 

.5624 

.6991 

.5850 

.6020 

.95 

.4895 

.3370 

.3510 

.6995 

.5805 

.5925 

.99 

.4895 

.3555 

.3429 

.6996 

.5773 

.5856 

Block  sizes  = ( 1,  5,  5,  3,  3,  2) 
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TABIZ  2 


2 

Variance  Reductions,  Rg  and  Rg,  for  Repairman  Problem 
Obtained  by  Simulating  the  Modified  Markov  Chain: 
r = E(X),  Return  State  = 0 


c 

Stopping  Time 

2 3 

s 

* 

^B 

1 

12 

.4676 

.3358 

.3602 

.6837 

.5795 

.6002 

2 

6 

.4597 

.3128 

.2321 

.6780 

.5595 

.5311 

3 

4 

.2950 

.2447 

.61^^ 

.5431 

4 

3 

.4476 

.2821 

.2432 

• 6691 

.5311 

•.»^952 

n = 10 
s = 4 
= 1 

* Block  sizes  = (1,  3,  5,  5,  3,  2) 


18 


4.  Cone lusions 

The  variance  reduction  technique  proposed  in  this  paper  shows 
promise  for  processes  with  relatively  long  cycle  lengths.  Examples  of 
such  processes  may  be  found  in  both  open  and  closed  queueing  networks. 

By  simulating  Che  modified  Markov  chain  (C^,  n > 0}  an  increased  number 
of  i.i.d.  cycles  are  obtained  for  a prespecified  run  length  chan  is 
possible  with  Che  original  chain  n > 0}.  Assuming  Chat  Che  chains 

n > 0)  and  n > 0)  are  approximately  equally  variable  over 

a cycle^  this  increase  in  the  number  of  cycles  translates  directly  into 
a variance  reduction.  Since  the  sample  paths  of  n > 0)  will 

usually  be  more  expensive  to  generate  than  those  of  (X^,  n > O),  the 
simulator  must  be  careful  to  ensure  that  Che  variance  reduction  obtained 
is  sufficient  to  produce  an  overall  computational  savings.  The  method 
has  the  best  chance  of  being  computationally  efficient  when  the  transition 
matrix  of  the  Markov  chain  is  relatively  sparse. 
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